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Pair cascades from millisecond pulsars (MSPs) may be a primary source of Galactic electrons and 
positrons that contribute to the increase in positron flux above 10 GeV as observed by PAMELA 
and AMS—02. The Fermi Large Area Telescope (LAT) has increased the number of detected y-ray 
MSPs tremendously. Light curve modelling furthermore favours abundant pair production in MSP 
magnetospheres, so that models of primary cosmic-ray positrons from pulsars should include 
the contribution from the larger numbers of MSPs and their potentially higher positron output 
per source. We model the contribution of Galactic MSPs to the terrestrial cosmic-ray electron / 
positron flux by using a population synthesis code to predict the source properties of present- 
day MSPs. We simulate pair spectra assuming an offset-dipole magnetic field which boosts pair 
creation rates. We also consider positrons and electrons that have additionally been accelerated to 
very high energies in the strong intrabinary shocks in black widow (BW) and redback (RB) binary 
systems. We transport these particles to Earth by calculating their diffusion and the radiative 
energy losses they suffer in the Galaxy using a model. Our model particle flux increases for non¬ 
zero offsets of the magnetic polar caps. We find that pair cascades from MSP magnetospheres 
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1. Introduction 


We now have firm evidence that the cosmic-ray positron fraction (PF) (j) (e+)/)], 
with (j) the flux, is an increasing function of energy above ~ 10 GeV, based on measurements by 
PAMELA Si, Fermi Large Area Telescope (LAT; [i), and the Alpha Magnetic Spectrometer 
{AMS—02-, iQ, |5[ 1^). Recent AMS—02 data extended the PF up to 500 GeV, indicating a levelling 
off of this fraction with energy, as well as being consistent with isotropy. If the PF is attributed 
solely to secondary positrons, produced during inelastic collisions between cosmic-ray nuclei and 
intergalactic hydrogen, this PF is expected to smoothly decrease with energy within the standard 
framework of cosmic-ray transport (e.g., 0)- The fact that the measured PF rises with energy 
may therefore point to nearby sources of primary positrons, either from dark matter annihilation 
(e.g., 0), or of astrophysical origin. The latter class of sources may include supernovae (e.g.,0), 
pulsar wind nebulae (e.g., [10]), young or mature pulsars (e.g., [|TT|, p^), and millisecond pulsars 


(MSPs; 0130). In this paper, we investigate the cosmic-ray flux contribution of the latter source 
class. 

MSPs are ancient pulsars that have been spun up to short rotational periods P by accretion from 
a binary companion QT^. The majority of MSPs were thought to lie below the pair creation death 
lines assuming dipole magnetic fields [|T5|], i.e., being ‘pair-starved’ [lAj. However, detection of 
narrow, double-peaked y-ray light curves trailing the radio peaks, very similar to those of younger 
pulsars, indicated the existence of narrow accelerator gaps in the MSP magnetospheres, requiring 
large numbers of electron-positron pairs to screen the electric field parallel to the magnetic field 
outside these gaps [|I7[]. Distortions of the surface magnetic field, either in the form of higher 
multipoles (e.g.. 


P) or offset polar caps (PCs; 119, 20]), may increase pair production in MSPs. 
MSPs furthermore produce electron-positron pairs with much higher energies than young pulsars 
due to their relatively low magnetic fields, so that the MSP pair spectra extend to several TeV [^]. 
Additionally, there has recently been a major increase in number of detected MSPs, many of them 
being nearby and relatively bright. MSPs are thus very promising potential sources of cosmic rays. 

In this paper we consider two MSP populations that may contribute to the terrestrial cosmic-ray 
flux. The first is a Galactic MSP population for which we obtain source properties via population 
synthesis modelling [^] (Section ^fl]). The second population involves MSPs with binary compan¬ 
ions. Shocks may form during collision of the the pulsar wind and the companion wind [ p^ p^ ], 
possibly accelerating the pairs escaping from the MSP magnetosphere to even higher energies. We 
therefore also consider black widow (BW) and redback (RB) systems as sources of cosmic rays 


(Section |2.2| ). The number of BWs and RBs has dramatically increased via Fermi observations. 
We calculate the source spectra originating in the MSP magnetospheres (Section |3.1[ ) and reac¬ 


celerated in binary shocks (for the BWs and RBs; Section 13.2] ). We next transport these spectra 
through the Galaxy to Earth (Section^ to assess the MSP contribution to the terrestrial cosmic-ray 
spectrum (Section Our conclusions follow in Section More details may be found in [^. 


2. Source populations 

2.1 Galactic synthesis model for the present-day MSP population 

). The Galaxy 


We follow [pi]] to predict the present-day distribution of MSPs (see also [ 25, 24 
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is seeded with MSPs (some aged up to 12 Gyr assuming a constant birth rate of 4.5 X 10^4 MSPs 
per century), which are then evolved in the Galactic potential from their birth location to the present 
time. We assume that MSPs are “born” on the spin-up line, and their surface magnetic field does 
not decay with time. We assume a power-law distribution for and adopt a spin-down power of 


the form (e.g., [26|) 




(2.1) 


where jJ. is the magnetic dipole moment, Q. is the rotational angular velocity, c is the speed of light, 
and a is the magnetic inclination angle relative to the pulsar’s rotational axis. We constrain the 
synthesis model parameters by fitting the model output to data from 12 radio surveys as well as 
from Fermi LAT. This simulation predicts the location as well as P and P (time derivative of P) of 
~50 000 Galactic MSPs, which we use as discrete sources of electrons and positrons. 

2.2 MSPs in binary systems - the BW / RB component 

About 80% of known MSPs are in binary systems. A subset of these, the BWs and RBs, 
may contain strong intrabinary shocks that can further accelerate the pairs. BWs are close binary 
systems with orbital periods of hours, containing a rotation-powered MSP and a compact com¬ 
panion having very low mass of ~ 0.01 — O.OSM©. The companion stars in BWs undergo intense 
heating of their atmospheres by the MSP wind, which drives a stellar wind and rapid mass loss 
from the star. A shock will form in the pulsar wind at the pressure balance point of the two winds 
and particle acceleration may occur in these shocks [^, ^]. RBs are similar systems, except that 
the companions have somewhat higher masses of ~ 0.1 — O.4M0 [27]. The MSPs in both types 


of system are typically energetic, with Lgd ~ 10^“^ — 10^^ergs Recent radio searches of Fermi 


unidentified y-ray point sources [28] have discovered many new BWs and RBs, with a total of 26 
known systems at the present time. By considering only the 24 publicly announced BWs and RBs 
here, our predictions will be a lower limit to the cosmic-ray flux contribution by binary MSPs. 


3. Source spectra 

3.1 Injection spectra of particles from the Galactic MSP population 

We calculate the spectra of pairs leaving the MSP magnetosphere using a code that follows the 


development of a PC electron-positron pair cascade in the magnetosphere [20[. A fraction of the 


curvature radiation photons emitted by primary particles ejected from the stellar surface undergo 
magnetic pair attenuation [p^. This produces a first-generation pair spectrum which then radiates 
synchrotron radiation (SR) photons that produce further generations of pairs. The total cascade 
multiplicity (average number of pairs spawned by each primary lepton) is a strong function of P 
and Bs- The sweepback of magnetic field lines near the light cylinder (where the corotation speed 
equals the speed of light) as well as asymmetric currents within the neutron star may cause the 


magnetic PCs to be offset from the dipole axis. We adopt a distorted magnetic field structure [2C] 
that leads to enhanced local electric fields, boosting pair formation even for pulsars below the usual 
pair death line. We consider PC offset parameter values e = (0,0.2,0.6) and use a grid in P and B^ 
to calculate the source spectrum for each source in the present-day MSP population (Section p4|) 
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via interpolation. From our simulations we find that about ~ 1 % of Lsd is tapped to generate the 
pairs. We neglect any further losses of the pair energy before injection into the interstellar medium, 
since MSPs are not surrounded by nebulae that can degrade the particle energy before escape. 

3.2 Injection spectra of particles accelerated in intrabinary shocks of BWs and RBs 

Pairs escaping from the pulsar magnetosphere may be further accelerated in the intrabinary 
shock that originates between the pulsar and companion winds in BW and RB systems. The maxi¬ 
mum particle energy will be determined by a balance between the minimum acceleration timescale, 
set by the particle diffusion (which we assume to be Bohm diffusion), and the SR loss timescale. 
This yields maximum particle energies in the TeV range [p^. We assume that the shock-accelerated 
spectrum will be an exponentially cut off power law with a spectral index of —2. We normalize 
this spectrum by requiring conservation of mass and energy (or equivalently, current and lumi¬ 
nosity). We assume a maximum shock efficiency T]p,max (conversion efficiency of Lsd to particle 
acceleration) of 10% and 30%. 


4. Galactic transport of injected leptons 


In order to transport the injected particles from the MSPs to Earth, we have to make some 
assumptions regarding the average Galactic background photon and magnetic field energy densi¬ 
ties. We approximate the interstellar radiation field using three blackbody components: optical, 
infrared (IR), and cosmic microwave background (CMB). We use two sets of energy densities, 
associated with the Galactic Disc and the Galactic Halo [p0|]. For the Disc, we assume f/opt = 
Gir = 0.4 eVcm^^ and Gcmb = 0.23 eVcm^^, while for the Halo we use Gopt = 0.8 eVcm^^, 
Gir = 0.05 eVcm^^, and Gcmb = 0.23 eV cm^^. For the average Galactic magnetic field strength 
that determines the SR loss rate we use values of Bsr ~ 1 — 3AtG ([Ull; hereafter DIO). We use a 
Fokker-Planck-type equation that includes spatial diffusion and energy losses: 


= V • (jr • VUe) - ^ (Ltotal^e) + S, 


(4.1) 


dt ^ dE 

with Ug the lepton density (per energy interval). Furthermore, denotes the diffusion tensor and 
Ltotai the total energy losses (we use SR and inverse Compton (IC) losses, the latter involving the 
full Klein-Nishina cross section), while S is the source term. Since MSPs are quite old (ages of 
~ 10^° yr), and have very small P values, we assume a steady-state scenario and invoke spherical 
symmetry. We assume that the diffusion coefficient is spatially independent so that becomes a 
scalar function of energy only: 

E 


k{E) = Ko 


ccd 


(4.2) 


We use typical values of Wd = 0.3 or 0.6, Lnorm = 1 GeV, and fCo = 0.1 kpc^Myr^' « 3 x 10^^ cm^s^ 
(e.g., 1^). For S, we consider At ~ 5 x 10“^ Galactic MSPs from the population synthesis code (Sec¬ 
tion |2.1| ), and N = 24 for the BW / RB case (Section |2!2| ). For the pulsar in our synthesis popu¬ 


lation, we assign a pair spectrum Qi{P,B^,e,E), as calculated in Section ^ for the corresponding 
simulated values of P, and £. We model this as 


N 


A = £a(^,fis,£)5(r-ro,,-). 


(4.3) 
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Figure 1: Total MSP contribution to the leptonic cosmic-ray spectrum at Earth, for Ko = 0.1 kpc^Myr^^ 
t?p,max = 0.1, Beff = 3.6 fJ.G, and f/opt = 0.4 eV cm^^. Electron spectra appear at the top and positron spectra 
at the bottom. The contribution from the synthesis component to the positron spectrum is (marginally) 
visible at ~ 30 GeV (for e = 0.6), and that of the BWs and RBs at ^ 1 TeV. Dashed and solid lines indicate 
curves for ajj = 0.3 and Uq = 0.6, respectively. The cool colours (purple, blue, and cyan) indicate spectra 
for e = 0.0,0.2, and 0.6. Green indicates the “background” (non-MSP) electrons and positrons predicted by 
GALPROP. 


Here, ro,, are the source positions. For a system of infinite extent, Equation ( |4.1[ ) is solved by the 
following Green’s function (e.g., lH]): 


G(r,ro,£',£'o) 


0(Eo-E) / 


r-ro|2\ 

A 


(4.4) 


with Eo the particle energy at the source, and the square of the propagation scale is given by 


A(E,Eo) 



<E’) 

£total(£') 


dE', 


(4.5) 


and @{Eq — E) the Heaviside function which ensures that A > 0. The lepton flux is given by 


U^,E) = G{r,ro,E,Eo)SdEod\o. 


(4.6) 


5. Results 

Figure [T] indicates the “background” secondary electron and positron fluxes predicted by GAL- 
PROP^ ||3^ for standard parameters, as well as data from PAMELA [^, Eertni ||3l], and AMS'—02 [|1]. 
We show our MSP synthesis spectral contribution (maximal but negligible contribution around 

* http://galprop.Stanford.edu/webrun/ 
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Figure 2: Measured ^ and predicted PF (including “background” contributions from GALPROP in 
green and the synthesis plus BW / RB contributions in purple, blue, and cyan, indicating e = (0.0,0.2,0.6). 
Here, tCo = 0.1 kpc^Myr^*, rjp max = 0.1, Beff = 3.6/rG, and t/opt = 0.4 eV cm^^. Dashed and solid lines 
are for = 0.3 and Ud = 0.6. 


30 GeV) plus BW / RB spectral contribution (we assume equal numbers of positrons and electrons) 
for dipole offsets of £ = (0.0,0.2,0.6), = (0.3,0.6), and (Beff,t^opt) = (3.6/iG,0.4 eVcm^^). 

Here, Beff includes the SR and Thomson-limit (non-optical) IC losses which are oc E'^. We set 
K() = 0.1 kpc^Myr^* and T]p max =0.1. The BW / RB contribution becomes higher for a larger 
T]p max> and all components are higher for lower values of Kq or (due to particle pile-up). The 
shape of the background model can strongly influence the total lepton spectrum. Figure shows 
the measured PF (e.g., |^) as well as the GALPROP and synthesis plus BW / RB contributions, for 
Ki) = 0.1 kpc^ Myr^^ and T]p,max = 0.1. The largest contribution is found ~ 100 GeV when e = 0.6 
and B = 3.6/iG. The highest PF (above 1 TeV) occurs for the lowest values of Kq and and 
highest flp^max- The BW / RB component makes a significant contribution at a few hundred GeV. 
Some parameter combinations are excluded by the data, e.g., Kq = 0.01 kpc^ Myr^^, rjp.max = 0-3, 
and £ = 0.6, depending on the background model. 

6. Conclusion 

We carefully assessed the contribution of MSPs to the cosmic-ray lepton spectra at Earth 
using a population synthesis code and a pair cascade code to calculate realistic source spectra. We 
also considered the contribution of binary BW / RB systems, which may further accelerate pairs 
escaping from the MSP magnetospheres in intrabinary shocks. The predicted MSP particle flux 
increases for non-zero magnetic field offset pai'ameters £. This is expected, since larger £ leads to 
an increase in the acceleration potential for some regions in magnetic azimuthal phase, implying 
an enhancement in both the number of particles as well as their maximum energy. The MSPs from 
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the synthesis model make only a modest contribution to the terrestrial cosmic-ray flux at a few 
tens of GeV, after which this spectral component cuts off. Increased magnetic and soft photon 
energy densities lead to increased particle energy losses in the Galaxy, and vice versa. The PF 
is somewhat enhanced above ~ 10 GeV by this component. The BW / RB component, however, 
contributes more substantially above several hundred GeV. For some parameter combinations, this 
component may even exceed the measured positron spectrum, and may violate the PF at high 
energies, depending on the background model. Alternative sources of primary positrons such as 
young, nearby pulsars or supernova remnants should also contribute to the cosmic-ray electron 
and positron flux. Future observations and modeling should continue to constrain the properties 
of these source classes, as well as improve our understanding of Galactic structure and particles 
within our Galaxy. 
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